#!/bin/bash

#PBS -l nodes=1:ppn=8,walltime=40:00:00
#PBS -N gatk_genotype
#PBS -M cmb433@nyu.edu
#PBS -m abe
#PBS -e localhost:/scratch/cmb433/awash_ddrad/rad-seq-primate-read-only/${PBS_JOBNAME}.e${PBS_JOBID}.${PBS_ARRAYID}
#PBS -o localhost:/scratch/cmb433/awash_ddrad/rad-seq-primate-read-only/${PBS_JOBNAME}.o${PBS_JOBID}.${PBS_ARRAYID}

# ------------------------------------------------------------------------------
# Variables
# ------------------------------------------------------------------------------

working_dir=/scratch/cmb433/awash_ddrad/rad-seq-primate-read-only

module load jdk/1.7.0

# ------------------------------------------------------------------------------
# Run pipeline
# ------------------------------------------------------------------------------

cd $working_dir

GATK=~/exome_macaque/bin/GATK
GENOME_FA=genomes/papAnu2/papAnu2.fa

if [ "$PBS_ARRAYID" -eq 21 ]; then
	CHROM=X
else
	CHROM=$PBS_ARRAYID
fi

BAMS=(`ls results/*.PE.bwa.baboon.passed.realn.bam`)

count=0
for b in ${BAMS[*]}; do
	BAMS[$count]="-I "$b" "
	count=`expr $count + 1`
done

java -jar ${GATK}/GenomeAnalysisTK.jar \
	-T UnifiedGenotyper \
	-R ${GENOME_FA} \
	${BAMS[*]} \
	-stand_call_conf 50.0 \
	-stand_emit_conf 10.0 \
	-o chr${CHROM}.raw.snps.indels.vcf \
	-nct 4 \
	-nt 8 \
	-L chr${CHROM}

exit;
